module history_mod

	use netcdf
	use namelist_mod
	use mesh_mod
	use state_mod
	use tracer_mod
	use vert_coord_mod, only : hyai, hybi

	implicit none

	private
	public history_init
	public history_output

contains 
	
	subroutine history_init()

		character(256) :: file_name
		integer ncid
		integer time_dimid, lon_dimid, lat_dimid, lev_dimid, ilon_dimid, ilat_dimid, ilev_dimid
		integer lon_varid, ilon_varid, lat_varid, ilat_varid, lev_varid, ilev_varid, cflx_varid, hyai_varid, hybi_varid
		integer q0_varid, q1_varid, q2_varid, q3_varid, u_varid, v_varid, we_varid, tm_varid, zs_varid
    integer q0_exact_varid, divx_varid, divy_varid
		integer cell_dim(4)

 	  character(len=*), parameter :: UNITS = "units"
 	  character(len=*), parameter :: lat_units = "degrees_north"
 	  character(len=*), parameter :: lon_units = "degrees_east"
 	  character(len=*), parameter :: lev_units = "downward"
 	  character(len=*), parameter :: ilat_units = "degrees_north"
 	  character(len=*), parameter :: ilon_units = "degrees_east"
 	  character(len=*), parameter :: ilev_units = "downward"

		file_name = trim(case_name) // '.nc'

		call check(NF90_CREATE(file_name, NF90_CLOBBER, ncid))

		call check(NF90_PUT_ATT(ncid, NF90_GLOBAL, 'scheme', 'FFSL'))

		call check(NF90_DEF_DIM(ncid, 'time', NF90_UNLIMITED, time_dimid))

		call check(NF90_DEF_DIM(ncid, 'lon', mesh%num_full_lon, lon_dimid))
		call check(NF90_DEF_VAR(ncid, 'lon', NF90_DOUBLE, [lon_dimid], lon_varid))
		call check(NF90_PUT_ATT(ncid, lon_varid, UNITS, lon_units))

		call check(NF90_DEF_DIM(ncid, 'ilon', mesh%num_half_lon, ilon_dimid))
		call check(NF90_DEF_VAR(ncid, 'ilon', NF90_DOUBLE, [ilon_dimid], ilon_varid))
		call check(NF90_PUT_ATT(ncid, ilon_varid, UNITS, ilon_units))

		call check(NF90_DEF_DIM(ncid, 'lat', mesh%num_full_lat, lat_dimid))
		call check(NF90_DEF_VAR(ncid, 'lat', NF90_DOUBLE, [lat_dimid], lat_varid))
		call check(NF90_PUT_ATT(ncid, lat_varid, UNITS, lat_units))

		call check(NF90_DEF_DIM(ncid, 'ilat', mesh%num_half_lat, ilat_dimid))
		call check(NF90_DEF_VAR(ncid, 'ilat', NF90_DOUBLE, [ilat_dimid], ilat_varid))
		call check(NF90_PUT_ATT(ncid, ilat_varid, UNITS, ilat_units))

		call check(NF90_DEF_DIM(ncid, 'lev', mesh%num_full_lev, lev_dimid))
		call check(NF90_DEF_VAR(ncid, 'lev', NF90_DOUBLE, [lev_dimid], lev_varid))
    call check(NF90_PUT_ATT(ncid, lev_varid, UNITS, lev_units))

    call check(NF90_DEF_DIM(ncid, 'ilev', mesh%num_half_lev, ilev_dimid))
		call check(NF90_DEF_VAR(ncid, 'ilev', NF90_DOUBLE, [ilev_dimid], ilev_varid))
		call check(NF90_PUT_ATT(ncid, ilev_varid, UNITS, ilev_units))

		call check(NF90_DEF_VAR(ncid, 'hyai', NF90_DOUBLE, [ilev_dimid], hyai_varid))
		call check(NF90_PUT_ATT(ncid, hyai_varid, UNITS, "hybrid coordinate coefficient"))
		call check(NF90_DEF_VAR(ncid, 'hybi', NF90_DOUBLE, [ilev_dimid], hybi_varid))
		call check(NF90_PUT_ATT(ncid, hybi_varid, UNITS, "hybrid coordinate coefficient"))


    cell_dim =[lon_dimid, lat_dimid, lev_dimid, time_dimid]

    call check(NF90_DEF_VAR(ncid, 'q0'      , NF90_DOUBLE, cell_dim, q0_varid))
    select case(test_case)
    case('solid_rotation')
      call check(NF90_DEF_VAR(ncid, 'q1'      , NF90_DOUBLE, cell_dim, q1_varid))
      call check(NF90_DEF_VAR(ncid, 'q2'      , NF90_DOUBLE, cell_dim, q2_varid))
    case('deform_case1','deform_case2','deform_case3','deform_case4')
	    call check(NF90_DEF_VAR(ncid, 'q1'      , NF90_DOUBLE, cell_dim, q1_varid))
	    call check(NF90_DEF_VAR(ncid, 'q2'      , NF90_DOUBLE, cell_dim, q2_varid))
	    call check(NF90_DEF_VAR(ncid, 'q3'      , NF90_DOUBLE, cell_dim, q3_varid))
	  case ('dcmip11')
	    call check(NF90_DEF_VAR(ncid, 'q1'      , NF90_DOUBLE, cell_dim, q1_varid))
	    call check(NF90_DEF_VAR(ncid, 'q2'      , NF90_DOUBLE, cell_dim, q2_varid))	
	  case('moving_vortex','stationary_vortex')
    	call check(NF90_DEF_VAR(ncid, 'q0_exact', NF90_DOUBLE, [lon_dimid, lat_dimid, time_dimid], q0_exact_varid))
    end select

    ! call check(NF90_DEF_VAR(ncid, 'cflx', NF90_DOUBLE, [ilon_dimid, lat_dimid , lev_dimid , time_dimid], cflx_varid))
    call check(NF90_DEF_VAR(ncid, 'u'   , NF90_DOUBLE, [ilon_dimid, lat_dimid , lev_dimid , time_dimid], u_varid   ))
    call check(NF90_PUT_ATT(ncid, u_varid, UNITS, "m s-1"))
    
    call check(NF90_DEF_VAR(ncid, 'v'   , NF90_DOUBLE, [lon_dimid , ilat_dimid, lev_dimid , time_dimid], v_varid   ))
    call check(NF90_PUT_ATT(ncid, v_varid, UNITS, "m s-1")) 
    if (mesh%num_full_lev > 1) then
    	call check(NF90_DEF_VAR(ncid, 'we'  , NF90_DOUBLE, [lon_dimid , lat_dimid , ilev_dimid, time_dimid], we_varid  ))
    	call check(NF90_PUT_ATT(ncid, we_varid, UNITS, ""))
    end if
    call check(NF90_DEF_VAR(ncid, 'tm' , NF90_DOUBLE, [time_dimid], tm_varid))
  
    call check(NF90_DEF_VAR(ncid, 'divx', NF90_DOUBLE, cell_dim, divx_varid))
    call check(NF90_DEF_VAR(ncid, 'divy', NF90_DOUBLE, cell_dim, divy_varid))

    if (test_case == 'dcmip13') then
    	call check(NF90_DEF_VAR(ncid, 'zs'      , NF90_DOUBLE, [lon_dimid, lat_dimid], zs_varid))
    	call check(NF90_PUT_ATT(ncid, zs_varid  , UNITS, "m"))
    end if

    call check(NF90_ENDDEF(ncid))

    call check(NF90_PUT_VAR(ncid, lon_varid , mesh%full_lon_deg(mesh%full_lon_ibeg:mesh%full_lon_iend)))
    call check(NF90_PUT_VAR(ncid, lat_varid , mesh%full_lat_deg(mesh%full_lat_ibeg:mesh%full_lat_iend)))
    call check(NF90_PUT_VAR(ncid, lev_varid , mesh%full_lev    (mesh%full_lev_ibeg:mesh%full_lev_iend)))
    call check(NF90_PUT_VAR(ncid, ilev_varid, mesh%half_lev    (mesh%half_lev_ibeg:mesh%half_lev_iend)))
    call check(NF90_PUT_VAR(ncid, ilon_varid, mesh%half_lon_deg(mesh%half_lon_ibeg:mesh%half_lon_iend)))
    call check(NF90_PUT_VAR(ncid, ilat_varid, mesh%half_lat_deg(mesh%half_lat_ibeg:mesh%half_lat_iend)))

    call check(NF90_PUT_VAR(ncid, hyai_varid, hyai))
    call check(NF90_PUT_VAR(ncid, hybi_varid, hybi))

    call check(NF90_CLOSE(ncid))

	end subroutine history_init

	subroutine history_output(state, tracer, time_step)

		type(state_type), intent(in) :: state
		type(tracer_type), intent(in) :: tracer
		integer, intent(in) :: time_step
		integer ncid
		integer q0_varid, q1_varid, q2_varid, q3_varid, u_varid, v_varid, we_varid, tm_varid, cflx_varid, zs_varid
    integer q0_exact_varid, divx_varid, divy_varid
		integer start(4), count(4)
		integer is, ie, js, je, ks, ke


		call check(NF90_OPEN(trim(case_name)//'.nc', NF90_WRITE, ncid))

		is = mesh%full_lon_ibeg; ie = mesh%full_lon_iend
		js = mesh%full_lat_ibeg; je = mesh%full_lat_iend
		ks = mesh%full_lev_ibeg; ke = mesh%full_lev_iend
		start = (/1,1,1,time_step/)
		count = (/mesh%num_full_lon, mesh%num_full_lat, mesh%num_full_lev, 1/)
    call check(NF90_INQ_VARID(ncid, 'q0', q0_varid))
		call check(NF90_PUT_VAR(ncid, q0_varid, tracer%q(is:ie,js:je,ks:ke,1), &
			         start=start, count=count))

		select case (test_case)
		case ('solid_rotation')
			call check(NF90_INQ_VARID(ncid, 'q1', q1_varid))
			call check(NF90_PUT_VAR(ncid, q1_varid, tracer%q(is:ie,js:je,ks:ke,2), &
			            start=start, count=count))
			
			call check(NF90_INQ_VARID(ncid, 'q2', q2_varid))
			call check(NF90_PUT_VAR(ncid, q2_varid, tracer%q(is:ie,js:je,ks:ke,3), &
			            start=start, count=count))
		case('deform_case1','deform_case2','deform_case3','deform_case4')
			call check(NF90_INQ_VARID(ncid, 'q1', q1_varid))
			call check(NF90_PUT_VAR(ncid, q1_varid, tracer%q(is:ie,js:je,ks:ke,2), &
								 start=start, count=count))

			call check(NF90_INQ_VARID(ncid, 'q2', q2_varid))
			call check(NF90_PUT_VAR(ncid, q2_varid, tracer%q(is:ie,js:je,ks:ke,3), &
								 start=start, count=count))
	    
	    call check(NF90_INQ_VARID(ncid, 'q3', q3_varid))
			call check(NF90_PUT_VAR(ncid, q3_varid, tracer%q(is:ie,js:je,ks:ke,4), &
				         start=start, count=count))	
		case ('dcmip11')
			call check(NF90_INQ_VARID(ncid, 'q1', q1_varid))
			call check(NF90_PUT_VAR(ncid, q1_varid, tracer%q(is:ie,js:je,ks:ke,2), &
								 start=start, count=count))

			call check(NF90_INQ_VARID(ncid, 'q2', q2_varid))
			call check(NF90_PUT_VAR(ncid, q2_varid, tracer%q(is:ie,js:je,ks:ke,3), &
								 start=start, count=count))
		case('moving_vortex','stationary_vortex')
    	call check(NF90_INQ_VARID(ncid, 'q0_exact', q0_exact_varid))
			call check(NF90_PUT_VAR(ncid, q0_exact_varid, tracer%qanl(is:ie,js:je), &
				         start=(/1,1,time_step/), count=(/mesh%num_full_lon, mesh%num_full_lat, 1/)))
		case default
		end select

		is = mesh%half_lon_ibeg; ie = mesh%half_lon_iend
		js = mesh%full_lat_ibeg; je = mesh%full_lat_iend
		ks = mesh%full_lev_ibeg; ke = mesh%full_lev_iend

		call check(NF90_INQ_VARID(ncid, 'u', u_varid))
		call check(NF90_PUT_VAR(ncid, u_varid, state%u(is:ie,js:je,ks:ke), &
		          start=(/1,1,1,time_step/), count=(/mesh%num_half_lon, mesh%num_full_lat, mesh%num_full_lev, 1/)))
		
		is = mesh%full_lon_ibeg; ie = mesh%full_lon_iend
		js = mesh%half_lat_ibeg; je = mesh%half_lat_iend
		ks = mesh%full_lev_ibeg; ke = mesh%full_lev_iend
		call check(NF90_INQ_VARID(ncid, 'v', v_varid))
		call check(NF90_PUT_VAR(ncid, v_varid, state%v(is:ie,js:je,ks:ke), &
		          start=(/1,1,1,time_step/), count=(/mesh%num_full_lon, mesh%num_half_lat, mesh%num_full_lev, 1/)))
	
		if (mesh%num_full_lev > 1) then
			is = mesh%full_lon_ibeg; ie = mesh%full_lon_iend
			js = mesh%full_lat_ibeg; je = mesh%full_lat_iend
			ks = mesh%half_lev_ibeg; ke = mesh%half_lev_iend	
			call check(NF90_INQ_VARID(ncid, 'we', we_varid))
			call check(NF90_PUT_VAR(ncid, we_varid, state%we(is:ie,js:je,ks:ke), &
			          start=(/1,1,1,time_step/), count=(/mesh%num_full_lon, mesh%num_full_lat, mesh%num_half_lev, 1/)))
	  end if

		call check(NF90_INQ_VARID(ncid, 'tm', tm_varid))
    call check(NF90_PUT_VAR(ncid, tm_varid, tracer%tm(1)))

		! is = mesh%half_lon_ibeg; ie = mesh%half_lon_iend
		! js = mesh%full_lat_ibeg; je = mesh%full_lat_iend
		! ks = mesh%full_lev_ibeg; ke = mesh%full_lev_iend
		! call check(NF90_INQ_VARID(ncid, 'cflx', cflx_varid))
		! call check(NF90_PUT_VAR(ncid, cflx_varid, tracer%cflx(is:ie,js:je,ks:ke), &
		          ! start=(/1,1,1,time_step/), count=(/mesh%num_half_lon, mesh%num_full_lat, mesh%num_full_lev, 1/)))
		
    if (test_case == 'dcmip13') then
    	is = mesh%full_lon_ibeg; ie = mesh%full_lon_iend
    	js = mesh%full_lat_ibeg; je = mesh%full_lat_iend
    	call check(NF90_INQ_VARID(ncid, 'zs', zs_varid))
    	call check(NF90_PUT_VAR(ncid, zs_varid, state%zs(is:ie,js:je), start=(/1,1/), count=(/mesh%num_full_lon, mesh%num_full_lat/)))
    end if

    is = mesh%full_lon_ibeg; ie = mesh%full_lon_iend
		js = mesh%full_lat_ibeg; je = mesh%full_lat_iend
		ks = mesh%full_lev_ibeg; ke = mesh%full_lev_iend
    call check(NF90_INQ_VARID(ncid, 'divx', divx_varid))
    call check(NF90_PUT_VAR(ncid, divx_varid, tracer%divx(is:ie,js:je,ks:ke), start=start, count=count))
    call check(NF90_INQ_VARID(ncid, 'divy', divy_varid))
    call check(NF90_PUT_VAR(ncid, divy_varid, tracer%divy(is:ie,js:je,ks:ke), start=start, count=count))

		call check(NF90_CLOSE(ncid))

	end subroutine history_output

	subroutine check(status)

		integer, intent(in) :: status

		if(status /= NF90_NOERR) then
			print*, trim(NF90_STRERROR(status))
			stop "Netcdf output stopped"
		end if

	end subroutine check

end module history_mod
